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ABSTRACT 

We explore the variability properties of long, high cadence GRMHD simulations across the elec¬ 
tromagnetic spectrum using an efficient, GPU-based radiative transfer algorithm. We focus on both 
disk- and jet-dominated simulations with parameters that successfully reproduce the time-averaged 
spectral properties of Sgr A* and the size of its image at 1.3 mm. We find that the disk-dominated 
models produce short timescale variability with amplitudes and power spectra that closely resemble 
those inferred observationally. In contrast, jet-dominated models generate only slow variability, at 
lower flux levels. Neither set of models show any X-ray flares, which most likely indicate that addi¬ 
tional physics, such as particle acceleration mechanisms, need to be incorporated into the GRMHD 
simulations to account for them. The disk-dominated models show strong, short-lived mm/IR flares, 
with short (< I hr) time lags between the mm and IR wavelengths, that arise from strong-field grav¬ 
itational leasing of magnetic flux tubes near the horizon. Such events provide a natural explanation 
for the observed IR flares with no X-ray counterparts. 

Subject headings: accretion, accretion disks — black hole physics — Galaxy: center — radiative 
transfer 


1. INTRODUCTION 

The black hole in the center of the Milky Way, Sgr A*, 
is the optimal laboratory for studying accretion physics 
at very low mass accretion rates. The mass and distance 
to the black hole are well constrained from dynamical 
studies of the orbits of stars in its v icinity (see, e.g., 
iGhez et ^ I2QQ8I: iGillessen et al.l l2QQ9f) and its bright¬ 
ness and variability have been studied extensively across 
the electromagnetic spe ctru m during the last decades 
(see iBaganoff et al.ll^QQII and iGenzel et al.ll2QQ3l for early 
studies). 

Long monitoring observations have revealed that 
Sgr A* is highly variable at timescales comparable to 
the dynamical time near i ts event horizon . In the in¬ 
frared (IR, e.g., K band), iDo et al.l (j2QQ9[ ) found that 
the power spectrum of variability can be modeled by 
a power law in the I-I5min range with an index that 
varies between —1.8 and —2.8. Similar power-law in¬ 
dices were obtained in other IR bands and at 230 GHz 

S er et al.ll2QQ8l:lYusef-Zadeh et ani2QQ9l : [Dexter et al.l 

. Moreover, iMever et al.l (|2QQ9[) found evidence for 
a break in the submillimeter ( submm) and IR p ower 
spectrum at ^ 150 min (see also lDexter et al.ll2QI4D . 

The amplitude of variability depends strongly on the 
wavelength of observation and changes from ^ 10% at 
1.3 mm to ^ 50% at 2/im, and u p to > 100% in the 
X-rays in three hour intervals (e.g.. iMarrone et ^120081: 
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iPoranet et al.l 1200^ iDo et al.ll2009[) . In the IR, it is de¬ 
bated whether the variability comprises of peaked fiar- 
ing events on top of a more quiescent variability (see 
iDodds-Eden et al]l2QIlL and references therein) or if all 
the observed variability is due to a continu um of events 
with a range of amplitudes and durations (jWitzel et al.l 
inil). Nevertheless, simultaneous IR/X-ray observa¬ 
tions of Sgr A* have established that even though all 
strong X-ra y flares have IR count erparts, the opposite 
is not true (|Hornstein et al.ll2QQ7l) . This suggests that 
either there are two physical origins for the IR events 
or that the physical mechanism responsible for the flares 
has two different observational manifestations. 

This wealth of observations has fueled substantial 
pro gress in modeling the a ccretion flow around Sgr A* 
fsee lYuan fc Naravanll2Qf3. for a review). In particular, 
the combination of general relativistic magnetohydro dy¬ 
namic (GRMHD) simulations and relativistic ray tracing 
calculations has led to a relative convergence of numer¬ 
ical models that account for numerous observed charac¬ 
teristics of Sgr A*, such as its quiesc ent broadband spec- 
trurn and the size of its mm image (iMoscibrodzka et al.l 
2QQ9I: [Scherbakov fc BaganoflDQIQl: iDexter et al.ll2QIQl: 
Chan et Sll2QI5[ ). However, a consensus has not been 
reached on the physical mechanisms responsible for the 
rich phenomenology of the variability observed from 
Sgr A* and, in particular, for the origin of the observed 
IR and X-ray flares. 

MHD simulations of accretion flows are in herently tur¬ 
bulent and variable but, as we showed in iChan et aTI 
(|2QQ9[) . the excursions in magnetic field and density gen¬ 
erated by MHD turbulence are not large enough to gen¬ 
erate the observed large amplitude o f variability , unles s 
additional physics is incorporated. In iChan et al l (j2QQ9f ). 
we investigated the effect of external perturb ations on 
the accretion flow. iDodds-Eden et alT (j2QIQ[ ) explored 
the possible connection between flares and the regions 
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of potential magnetic re connection in MHD simulations. 
iDexter fc Fragilel (j2Q13f ) investigated the possibility of 
shock heating in a tilted accretion disk as an alternate 
avenue to generating large variability and flares. 

In parallel to the numerical simulations, a num¬ 
ber of phenomenological models of the variability have 
been explored, which involve, e.g., orbiting hotspots 
or sudden inject ion of non-thermal particles (e.g., 
Markoff et al.l[2QQ^ lYnan et al.ll2QQ4 iDodds-Eden et al.l 

2QQ9I : lEckart et al.ll2QQ9l: iDibi et al.ll2Q14l) . These models 
can reproduce some characteristics of the data but their 
connection to the turbulence and thermodynamics of the 
accretion flow is unclear. 

The developments in the algorithms and the compu¬ 
tational power in the last two years allow us to revisit 
the question of the broadband Sgr A* variability us¬ 
ing long GRMHD simulations, with high cadence ray 
tracing to compute observables, and more detailed de¬ 
scriptions of the thermodynamics of the accretion flow. 
We make use of GRMHD simulations o f the accretion 
flow around Sgr A* perform ed with HARM (jNaravan et "ahl 
l2Q12l : iSadowski et al.l[2QT^ and radiative transfer calcu¬ 
lations performed with the highly efficient G PU based 
ray-tracing algorithm GRay (jChan et al.ll2Q13l) . In earlier 
work, we used the combination of these two algorithms 
to carry out a comprehensive study of models with a 
wide range of black hole spins, magnetic field geome¬ 
tries, plasma physics models, and geometric parameters 
at hig h time, energy, and spatial resolutions (jGhan et al.l 
120151 ). Among these, we identified five simulations with 
different characteristics that all reproduce the observed 
quiescent broadband spectrum and the 1.3 mm image size 
of Sgr A*. We use these as our five baseline numerical 
models. 

In this paper, we calculate multi-wavelength 
lightcurves with a 10 GMc“^ (« 3.5 min for Sgr A*) 
time resolution that span approximately 10, 000 GMc“^ 
(60 hr) for each of the five simulations. We use these 
lightcurves to compute the amplitude of variability 
and the power spectrum over a wide range of wave¬ 
lengths and to investigate the correlations and potential 
lags between different wavelengths. We find that the 
disk-dominated (SANE) models reproduce the short 
timescale mm and IR variability with amplitudes and 
power spectra that are comparable to those inferred 
observationally, whereas jet-dominated (MAD) models 
produce significantly less variability. We also find that 
none of the models produce X-ray flares, which most 
likely require additional physics, such as particle accel¬ 
eration, that is not captured within the simulations. We 
trace the origins of the large amplitude IR variability 
to gravitational leasing of short-lived, bright flux tubes 
near the black hole horizon, which produces achromatic 
peaks in the lightcurves. We discuss the prospect of 
identifying these leasing events with future imaging 
experiments at mm and IR wavelengths. 

2. THE GRMHD+RAY TRACING SIMULATIONS 

In I Chan et al.l (j2Q15[ ) we considered a large number of 
GRMHD models, with a wide range of black-hole spins 
(0 to 0.9), accretion rates (which determine the scale of 
electron density in the accretion flow), initial configu¬ 
rations of the magnetic field (SAN E and MAD in the 
terminology of iNaravan et al.ll2Q12f ). and with different 
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Figure 1. The mean spectra of the five best-fit GRMHD models 
to Sgr A* (see Table m and their observational constraints. The 
colorful curves are the mean spectra iy{Lu) computed using all 1024 
snapshots. The purple curve is plotted in dash for clarity when it 
overlaps with the blue or green curves. The crosses, vertical ar- 
ro ws, and bow-ties are the observational data that we employed 
in ICtian et al.l (120151 1 to constrain the models. The gray band be¬ 
tween u ^ 10^^ Hz and 10^^ Hz marks the frequency range over 
which we perform a least-squares fit to the radio spectrum. The 
gray line ai u ^ 10^^ Hz marks the IR frequency at which we used 
the range of fiuxes observed at different times to impose an upper 
and lower bound on the models. The gray line at u ^ 10^^ Hz 
marks the X-ray frequency, at which we used 10% of the observed 
quiescent ffux (the open circle below the bowtie) to fix the density 
normalization in the fiow. The white dotted line inside the gray 
band marks A = 1.3 mm, at which we fit the models to the EHT 
measurement of the image size. 

prescriptions for the thermodynamic properties of the 
electrons in the accretion flow and in the magnetically 
dominated funnel. We then used GRay to compute av¬ 
erage images and spectra from 100 snapshots from each 
simulation, with a time resolution of lOGM/c^, where 
G is the gravitational constant, M = 4.3 x 10^ Mq is 
the mass of the Galactic black hole, and c is the speed 
of light. Hereafter, we will use gravitational units with 
G = c=l. 

We identified five GRMHD models that are able to ac¬ 
count for the average spectrum of Sgr A* and the mea¬ 
sured size of its emission region at 1.3 mm, as measured 
by the EHT. In particular, we required that each success¬ 
ful model reproduces (a) the correct average flux and 
spectral slope at lO^^-lO^^Hz, (b) a. flux at 10^^ Hz 
that falls within the observed range of highly variable 
IR fluxes, (c) an X-ray flux equal to 10% of the observed 
quiescent flux from Sgr A*, which has been attributed by 
iNeilsen et al.l (j2Q13l ) to emission from the inner accretion 
flow, and (d) a scatter-broadened image at 1.3 mm with 
a size th at is equal to the charac teristic size measured by 
the EHT iDoeleman et al.l (|2QQ8f ). We list the parameters 
of these five models in Table [H 

In this paper, in order to study the variability predicted 
by the GRMHD simulations, we use 1024 late-time snap- 
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Table 1 _ 

Parameters Of The Five GRMHD Models from I Chan et all (120151) 


Model 

a 

Bo 

^e(^disk) 

^disk 

^e,funnel 

^funnel 

^begin 

^end 

A 

0.7 

SANE 

6.9 X lO'^ 

0.02371 

56.2 

— 

89,760 

99,990 

B 

0.9 

SANE 

5.5 X 10'^ 

0.01000 

31.6 

— 

39,760 

49,990 

G 

0.0 

MAD 

5.9 X 10^ 

0.00056 

— 

0.0100 

209,760 

219,990 

D 

0.9 

MAD 

1.6 X 10® 

0.00075 

1.8 

— 

38,520 

48,750 

E 

0.9 

MAD 

1.6 X 10® 

0.00075 

— 

0.0032 

38,520 

48,750 


Note. — The parameters of the five best-fit GRMHD models of Sgr A* that we obtained in IChan et al.l (120151 ). 
The first column lists the model names that we use throughout this paper. The parameters a, Hq, Ue, and ^disk are 
the black hole spin, the initial magnetic field configuration, the normalization of the electron number density, and 
the electron-ion temperature ratio in the accretion fiow. The parameters Te,funnel and ^funnel are the dimensionless 
electron temperature and electron-ion temperature ratio in the funnel regions. The times listed in the last two 
columns correspond to the initial and final simulation times, in gravitational units, GMc~^, that we use for the 
study of variability presented here. Note that there are 1024 snapshots for each model separated by lOM, so that 
^end “ ^begin is always equal to 10230GMc“^, which is approximately 60hr for the mass of Sgr A*. 


shots of the best-fit models, with the same time resolu¬ 
tion (lOM). For the mass of Sgr A*, this corresponds 
to a time resolution of ~ 212 s with a span of ~ 60 hr. 
For each snapshot, we compute four levels of images at 
each wavelength, with sizes equal to 1024M x 1024M, 
256M X 256M, 64M x 64M, and 16M x 16M. Each 
level contains 512 x 512 pixels. Hence, the resolution im¬ 
proves by a factor of 4 in each direction going from one 
level to the next. We then calculate the total flux at each 
time and at each wavelength by composing together all 
four levels of images. In presenting the resulting flux, 
we often use the quantity which would 

have been equal to the luminosity of the source, had the 
emission been isotropic. 

In the ray tracing calculations, we use the fast 
light approximation^ where the speed of light is as¬ 
sumed t o be infinite, as was done in numerous previous 
studies (iDolence et al.l I2QQ9I: iMoscibrodzka et I2Q11I: 
iShiokawa et al.|[2Q12[) . This approximation is valid when 
the light crossing time ^crossing = Rjc across a region of 
size R is much less than the dynamical time scale at that 
radius tdynamicai = i.c., when R > GMjc?. 

This is true for the emission from Sgr A* for wavelengths 
longer than ^ mm or shorter than ^ /am (IR). At mm 
to IR wavelengths, where the emitting region has a size 
that is comparable to the horizon of Sgr A*, the light 
crossing time is ~ 20 s, which is a factor of ten smaller 
than the time resolution of the snapshots. 

In Figure [TJ we plot the average predicted spectrum for 
each m odel together with the observational constraints 
used in iChan et "aTI (|2Q15[) . Even though we have ex¬ 
panded the time span of the snapshots used by a factor 
of 10, the average spectra do not violate any of the ob¬ 
servational constraints. 

3. THE VARIABILITY OF SIMULATED LIGHTGURVES 

In this section, we study the lightcurves from each of 
the five models using a variety of time-series analysis 
techniques. We compute the power spectra, the frac¬ 
tional rms amplitudes as a function of photon frequency 
and investigate correlations between different frequen¬ 
cies. 

3.1. Lightcurves 

We first compute the lightcurves Ly(t) from each of 
the five models at four representative frequencies. In each 


panel of EigureO we plot the radio {u = 10^^ Hz), 1.3 mm 
(z/ = 2.31 X 10^^ Hz), IR (z/ = 1.38 x 10^"^ Hz), and X-ray 
[y = 1.04 X 10^^ Hz) lightcurves in red, black, green, and 
blue, respectively. We will use this color scheme when 
we compare different wavelengths throughout the paper. 

The fluxes and the variability amplitudes span a wide 
dynamical range across the wavelengths of interest, as 
we will discuss below. In order to be able to compare 
the lightcurves at different wavelengths, we first subtract 
the mean flux at each wavelength from the correspond¬ 
ing lightcurve and normalize it to the standard deviation 
calculated over the entire interval of the simulation. As 
a result, the mean is zero and the standard deviation 
is equal to unity for the normalized flux at each wave¬ 
length. In the X-rays, there is also an overall reduction 
of flux with time because of the depletion of the torus 
that is feeding the black hole. (This is pronounced in 
the X-rays because the emission originates from a large 
volume). To remove this artifact of the simulation setup, 
we also detrend the X-ray lightcurve by fitting a line to 
the X-ray flux as a function of time and subtracting this 
from the lightcurve. 

We show in Eigure[3]the fractional rms amplitudes over 
the entire length of the simulations as a function of pho¬ 
ton frequency for the five different models, which are 
the standard deviations that we used to normalize the 
lightcurves shown in EigureO At low photon frequencies, 
the amplitudes increase with frequency, peaking near the 
optical at ~ 10^^ Hz. This is the transition between the 
synchrotron and the bremsstrahlung dominated parts of 
the spectrum. At the wavelengths where the synchrotron 
emission dominates, a small change in the density, tem¬ 
perature, or the magnetic field can cause a large change 
in the flux and, therefore, give rise to a large amplitude 
of variability. In contrast, the Bremsstrahlung emission 
is optically thin and originates from large volumes, aver¬ 
aging out fluctuations in density and temperature. 

We can draw several conclusions from Eigures [2] and 
[3l Eirst, the jet-dominated models, which are the MAD 
models D and E, overall have larger amplitudes of vari¬ 
ability but their radio-to-IR lightcurves are dominated 
by smoother, long timescale flux changes. In contrast, 
the disk dominated models (SANE Models A and B), 
have overall smaller amplitudes of variability but their 
lightcurves are characterized by short-lived, high ampli¬ 
tude excursions in their IR and mm flux. As we will 
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Figure 2. Lightcurves of the best-fit models. In each panel, we plot the radio, 1.3 mm, IR, and X-ray lightcurves in red, black, green, and 
blue, respectively. We offset and normalize all lightcurves for easy comparison. In addition, we also detrend the X-ray lightcurves. 


discuss in more detail below, these excursions have a 
lot of the characteristics of the IR flares observed from 
Sgr A*. Model C, in which the long-wavelength emis¬ 
sion originates primarily from the base of the jet, has 
characteristics of both disk- and jet-dominated models 
and shows both long and short timescale variability. At 
1.3 mm, the Model C lightcurve shows evidence for quasi- 
periodic dimming, the origin of which we will explore be¬ 
low. Finally, in the X-rays, all models show very weak 
and smooth variations in the flux, with no evidence for 
any flaring events such as those observed from Sgr A* by 
Chandra. 

Our multi-wavelength lightcurves allow us to look for 
correlations between the variability at different wave¬ 
lengths. In the jet-dominated models, all wavelengths 
(except X-rays) show quite similar trends in variability 
and are roughly synchronized. There is some weak evi¬ 


dence of the radio lightcurve lagging behind the 1.3 mm 
and the IR (see e.g., the times between 6000 M and 
8000 M in the lower panels of Figure [2]). In the disk- 
dominated models, it is only the 1.3 mm and the IR 
lightcurves that vary in sync, while the radio is decou¬ 
pled. Moreover, the sharp IR flares have nearly coinci¬ 
dent, albeit weaker, 1.3 mm counterparts with only very 
short (< 1 hr) time lags. In all models, the X-ray vari¬ 
ability is decoupled from the long wavelengths, which is 
expected given their very different emission locations in 
the flow. 

3.2. Flux Histograms 

The range and character of variability observed from 
Sgr A* is often investigated via studies of flux histograms 
(see, e.g.. lYusef-Zadeh et al.ll2QQ6l. I2QQ9I: iDo et al.l]2QQ9l: 

iDodds-Eden et aH l2Qlll: IWitzel et al.l l2Ql^. Figure 
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Figure 3. The fractional rms amplitude computed using the full 
lightcurve. The color scheme of the amplitudes is identical to the 
one used in Figure[T] except that the purple curve is plotted in dash 
when it overlaps with the blue curves. The amplitudes peak near 
the visible frequencies, i.e., u ~ 10^^ Hz, which is the transition 
region between the synchrotron- and bremsstrahlung-dominated 
regions. 

shows the flux histograms at 2.2 fim for our five mod¬ 
els over the entire simulation time, with each flux corre¬ 
sponding to a time integration of lOM. (We chose this 
wavelength for thi s plot, in order to compy e it directly 
to the analysis of iDodds-Eden et al.ll2Qlll ) As in the 
case of the observational studies, we fit the low-flux ends 
of the histograms (up to 0.6 x lO^^ergs”^ for Models A 
and B and up to 0.03 x 10^^ ergs“^ for Models C-E) with 
a log-normal distribution. 

The flux histograms of the disk-dominated SANE Mod¬ 
els A and B are well described by the log-normal distri¬ 
bution with marginal evidence for a small excess at large 
fluxes. In other words, the same physical mechanism ap¬ 
pears to be responsible for most of the variability, with 
some evidence that the brightest of flares occur at larger 
than expected fluxes. As we will discuss in the next 
section, this is consistent with our interpretation that 
most of the variability in these models is the result of 
the presence of short-lived hot, magnetically dominated 
flux tubes in the accretion flow. However, the brightest 
flares are generated when one of these magnetically dom¬ 
inated regions gets lensed behind the black hole and its 
brightness is greatly magnified. 

The flux histograms of the jet-dominate MAD Mod¬ 
els C-E are instead significantly flatter than the log¬ 
normal distribution at all but the smallest fluxes. This is 
another illustration of the fact that the IR lightcurves of 
these models show slow variability over a range of fluxes 
that cannot be described in terms of a range of short¬ 
lived flares. 

3.3. Fractional RMS Amplitudes 


In Eigure [3l we looked at the frequency dependence 
of the fractional rms amplitude using the entire ^ 60 hr 
time span of the simulations. However, given the contin¬ 
uum of flux excursions and timescales visible in Eigure O 
these amplitudes will depend on the time span used to 
calculate the rms and the particular interval for which it 
is computed. This is especially important when compar¬ 
ing theoretical amplitudes to observations, which usually 
span intervals that are much shorter. 

In order to carry out the comparison to observations, 
we use th e simultaneous multi-w avelength monitoring of 
Sgr A* of iMarrone et al.l (j2QQ8[ ). which is one of the few 
data sets that spans the wavelengths from radio to X- 
rays. In particular, we use the 2006 July 17 lightcurve, 
during which an X-ray flare was detected. The longest 
time span of simultaneous observations is equal to 3 hr 
and the range of rms amplitudes measured over a num¬ 
ber of overlapping 3 hr intervals during this observation 
is shown as red data points with error bars in Eigure [5l 
In the same figure, we plot as gray lines the predicted 
fractional rms amplitudes computed from many overlap¬ 
ping 510 M (^ 3 hr) intervals from the simulations. 

As expected, the simulations show a wide range of frac¬ 
tional amplitudes at each wavelength depending on the 
particular realization of the turbulent magnetic field and 
density in the accretion flow. However, all realizations 
have the same overall wavelength dependence and peak 
at optical frequencies. Eurthermore, the data points in 
the mm and IR overlap with typical realizations of the 
simulated curves. On the other hand, the X-ray variabil¬ 
ity amplitude, which reflects the presence of the X-ray 
flare, is many orders of magnitudes larger than any of 
the predicted realizations. 

This comparison between the observed and simulated 
variability amplitudes indicates that GRMHD simula¬ 
tions naturally produce IR flares and variability with no 
X-ray counterparts. Indeed, observations show two types 
of flaring events; IR flares are often not accompanied 
by X-ray counte rparts, whereas X-ray flares always have 
IR counterparts (|Hornstein et al.l[2QQ7[ ) . Interpreting the 
observations in the context of our simulations suggests 
that there are two distinct mechanisms that cause the 
two types of flaring activity. As we will discuss in detail 
in the next section, strong-field gravitational leasing in 
the turbulent GRMHD flows naturally produces the ma¬ 
jority of IR flares without X-ray counterparts. However, 
additional physical processes must be incorporated into 
the simulations to generate the X-ray flares. 

Given the sensitivity of the variability amplitude on 
the time span used, we can use the simulations to ex¬ 
plore this dependence in detail. In Eigure [6l we show 
the fractional rms amplitude at two different mm and IR 
wavelengths each, as a function of the time interval over 
which it is evaluated. Eor each value of the time inter¬ 
val, the red points correspond to the average fractional 
rms amplitude for the same overlapping segments in the 
simulation discussed earlier and the error bars represent 
the variance among them. As expected, the overall trend 
is that the amplitude increases with increasing time in¬ 
terval, with no evidence for a characteristic time scale in 
the range explored (i.e., 5-40hr). This is typical of red 
noise processes, which we can explore in more detail us¬ 
ing power spectral techniques. In addition, the fact that 
the fractional rms amplitude in the IR becomes larger 
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Figure 4. IR flux histograms for the flve models, evaluated at 2.2/xm (the green curves in Figure[2j, together with the best-fit log-normal 
distributions. The disk-dominated SANE Models A and B are well described by the log-normal distribution, with marginal evidence for a 
small excess at large fluxes. Noting the different scales in the horizontal axes, the jet-dominated MAD Models C-E span a much narrower 
range in fluxes but their distribution is significantly flatter than the log-normal curve. 


than 100 % points to a large number of short-lived, high 
amplitude flux excursions (flares) being responsible for 
this variability. This observation also explains the fact 
that the range of variability amplitudes in different 3 hr 
segments is larger in the IR than at mm wavelengths, 
since the amplitude of any particular segment in the IR 
depends strongly on whether a strong flare has occurred 
within its duration. 

3.4. Power Spectra 

Red noise variability is best studied in Fourier space, 
using power-spectral techniques. However, generating a 
power spectrum requires uniformly sampled light curves, 
with fast sampling and minimal data gaps. This is very 
difficult to achieve with ground-based observations of 
Sgr A* (e.g., in the mm and in the IR), which would 
require data spanning many hours but sampled every 
few minutes. For this reason, very few attempts have 
been made to measure the power spectrum of Sgr A*, 
and the few efforts have primarily been in the IR via the 
measurement of the structure function. In all cases, the 
power-spectrum was modeled by a power-law function 
(iMeyer et al.l l2QQ8l iDo et al.l 120091 IZamaninasab et 
2Q1QL andlYusef-Zadeh et"am2QQ9l for near-infrared, and 
Dexter et al.ll2Q14l for submmb In a very detailed study 
of the systematic uncertainties that are inherent to the 
measurement o f power spectra via structure functions, 
iDo et all (j2QQ9f ) reported indices for the K band power 
spectr a in the range —1.8 to —2.8, while iMever et al.l 
(j 2 QQ 8 f ) reported evidence for a break in the power 
spectrum at minutes, with the power-law index 

changing from —2.1 ± 0.5 above this break to — 0 . 3 lo '4 
below it (see also lDexter et al.1l2Q14l for a similar report 
of a break at 8^4 hr for submm wavelengths). 

To make a quantitative comparison of the variability 
predicted in the GRMHD simulations to the variability 
reported in these studies, we compute the power spec¬ 
tra for each of the five models at four different photon 
frequencies and plot them in Figure [71 To compute the 
power spectra, we divide the 1024 snapshots from each 
light curve into 7 overlapping segments, each of which 


contains 256 data points. We multiply these segments 
by a symmetric “triangle” window function. We then 
compute the fast Fourier transform of the windowed seg¬ 
ments and average them. Finally, we bin the spectra with 
a bin size equal to the Fibonacci sequence. Each panel of 
Figure [7] contains the power spectra for the five models 
A-E, color-coded as before, for the radio, 1.3 mm, IR, 
and the X-ray light curves. 

The power spectra in the X-rays (fourth panel in Eig- 
ure El) show an / ^ dependence, which is characteris¬ 
tic of long-wavelength variability over timescales that 
are longer than the span of the simulations. We can, 
therefore, use them only to infer that there is no evi¬ 
dence for fast X-ray variability in the simulations. In 
all three of the longer wavelengths, however, the models 
give rise to power spectra with a range of flatter slopes. 
In the radio, the slopes can be approximately charac¬ 
terized as with Models C and E having somewhat 
steeper power spectra. At 1.3 mm, all models show a 
flatter dependence, with slopes comparable to « — 2. In 
addition. Model C has the flattest spectrum among all 
of the models, consistent with the fact that it has the 
strongest short-timescale variability at this wavelength 
(see Eigure|2]). 

In the IR, the two disk-dominated models (A and B) 
produce higher relative power at higher frequencies, as 
expected from the visual inspection of lightcurves that 
show high levels of short timescale variability. These two 
models also exhibit some weak evidence for flattening to¬ 
ward smaller frequencies. To compare these power spec¬ 
tra to observations, we also plot in dashed line a power 
spectrum that has a break at 154 min, with slopes —0.3 
and — 2.1 below and above the break and a gray band 
that corresponds to t he formal uncertai nties in these 
slopes (as inferred by IMever et "ahl I2QQ9D . This figure 
demonstrates that not only the overall amplitude but also 
the detailed power spectrum of variability generated by 
our disk-dominated GRMHD simulations are consistent 
with Sgr A* observations. Therefore, the observed IR 
power spectra appear to favor the disk-dominated SANE 
models as opposed to the jet-dominated MAD ones. 
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Figure 5. Gray curves show the fractional rms amplitudes computed usiug uumerous overlappiug 510M 3 hr) iutervals from the 

simulatious, startiug at times t — tbegin = 0? lOOM, 200M,.... The blue curve iu each pauel, which shows the rms amplitude as a fuuctiou 
of photou frequeucy computed usiug the eutire spau of the simulatiou, is ideutical to the correspoudiug curve iu Figure [3] The red data 
poiuts aud error bars m ark the rauge of rms amplitudes measured over a uumber of overlappiug 3 hr iutervals duriug the 2006 July 17 
observatious reported iu lMarroue et all (|2008l) . 


3.5. Correlations 

A higher order Fourier statistic that is often used in 
the variability studies of Sgr A* is the cross correlation 
function between different wavelengths that is used to 
search for time lags in the lightcurves. Such lags are 
often interpreted as signatures of short-lived, hot emit¬ 
ted regions that are advected with the accretion flow or 
are ejected in a je t or outflow (e.g. Mar rone et al.l I2QQ8I: 
lEckart et al.ll2Q()9l : iBrinkerink et al. I2Q15D . We compute 
the cross correlation between the 1.3 mm lightcurve and 
the lightcurves at other wavelengths to look for such sig¬ 
natures. 

The first three panels in Figure [HI from left to right, 
show the cross correlations as functions of time lag be¬ 
tween the 1.3 mm lightcurve and the radio, IR, and X-ray 
lightcurves, respectively. A positive lag implies that the 
lightcurve at the particular wavelength lags behind the 
1.3 mm lightcurve, whereas a negative lag implies that it 
leads the 1.3 mm lightcurve. There are clear differences 
between the 1.3mm-radio correlations and 1.3mm-IR 
correlations with the 1.3mm-X-ray correlations. This is 
not surprising, because most of the X-ray emission comes 
from large radii, while the 1.3 mm emission originates 
near the event horizon. In order to focus on the 1.3 mm- 
X-ray correlations of the inner flow, we also compute X- 
ray lightcurves using the inner 16M x 16M images only, 
and cross correlate the 1.3 mm lightcurves with these in¬ 
ner X-ray lightcurves. The result is plotted in the fourth 
panel of Figure |8l 

In the disk-dominated, SANE, models, there is clear 
evidence for correlated variability only between 1.3 mm 
and IR curves with small (< 1 hr) time lags. In Sec¬ 
tion 14.1.11 we use the simulated images to argue that 
variability due to strong-held gravitational leasing is the 
cause for the achromatic character of the large hares seen 
in the simulations as well as for the lack of time lags. On 
the other hand, there is no correlated variability between 
radio, 1.3 mm, and X-ray wavelengths, as can also be seen 
in the lightcurves shown in Eigure|2l The two peaks ap¬ 
pearing at time lags ^ 10 hr for models A and B between 


the radio and 1.3 mm as well as the two peaks at ^ —10 hr 
time lags between the inner X-rays and the 1.3 mm sim¬ 
ply correspond to the time diherences between physically 
unrelated peaks in the lightcurves. Such chance events 
emphasizes the difficulty of using cross-correlation anal¬ 
ysis in short observations to infer the time evolution of 
emission structures in the accretion flows. 

In the jet-dominated, MAD, solutions, the radio, 
1.3 mm, and IR lightcurves are highly correlated, with 
marginal evidence of radio lightcurves lagging and the 
IR leading the 1.3 mm lightcurves by ^ 1 hr. There also 
seems to be a similar level of correlation between the in¬ 
ner X-ray and 1.3 mm lightcurves, with the former lead¬ 
ing the latter by ^ 1 hr. This ordering suggests that the 
physical mechanism driving the variability in these so¬ 
lutions causes first the X-rays (in the inner disk), then 
the IR, then the 1.3 mm, then finally the radio fiux to 
increase. In section H. 2. 2[ by studying the simulated im¬ 
ages, we conjecture that magnetic reconnection is respon¬ 
sible for the variability in these models. 

4. THE PHYSICAL ORIGINS OF THE VARIABILITY 

In the previous section, we focused on calculating 
lightcurves of the simulated emission at different wave¬ 
lengths and their statistical properties. In this section, 
we take a deeper look at the simulated time-dependent 
images in order to explore and understand the physical 
origins of the different types of variability that we see in 
the simulations. 

4.1. Disk-dominated SANE models 

The variability in the disk-dominated SANE models is 
characterized by a spectrum of intense, short lived flares 
in the IR, with related (albeit weaker) flux excursions 
at 1.3 mm, and uncorrelated slow variability in the radio 
and in the X-rays. 

We first discuss the origin of the IR/1.3mm flares, 
which we attribute to strong-held gravitational leasing 
of hot hux tubes in the accretion how and then explore 
the origin of the radio variability, which we attribute to 































10^ 


10^ r 


t 


}}n»‘ 


i • 


1.3mm 
Model A 


H-^^^^-h 


} i M H ‘ • 


i n 


f • 


0.84mm . 
Model A 

^^- 


IIIIIH 


i • 


4.5f^m 

Model A 

H-^^- 


jtlfiH 


i * 


2.1/im 
Model A 
H-^^- 


10 r 


P 10^ 






1.3mm 
Model B 


H- ^^^^-h 


H H ‘« ‘uL' 

I Model C 


}}HH» 


0.84mm 
Model B 


IfH 


i i 


4.5pim 

Model B 


flti 


2.1^m 
Model B 


I 10 

cd 


cd . ^1 
10 
o 


^ 2 
10 


I i 0.84mm 


H- ^^^^-h 


Model C 
H-^- 


H-h 




4.5;fm 

Model C 


M 


}• 


2.1/im 
Model C 


H-^-K 


10' 


• I 1.3mm 

Model D 

-H - ^^^^^^- 


1 I 1.3mm 

I Model E 


f 1 0.84mm. 

1 Model D 


}}’■ 


4.5f4.m 

Model D 


M 


2.1^m 
Model D 


10^ r 


10' r 


ilHfn- 

1 ^ 0.84mm 

I Model E 


}}»■ 


4.5;fm 

Model E 


{{ 


2.1/im 
Model E 


0 5 10 15 20 25 30 35 0 5 10 15 20 25 30 35 0 5 10 15 20 25 30 35 0 5 10 15 20 25 30 35 

Time interval (hr) 


Figure 6. The fractional rms amplitude at 1.3 mm, 0.86 mm, 4.5/xm, and 2.1/xm, as a function of the time interval over which it is 
evaluated. For each value of the time interval, the red points correspond to the average fractional rms amplitude for all the overlapping 
segments in the simulation starting at times t — tbegin = 0? lOOM, 200M,. .. and the error bars represent the variance among them. 
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Figure 7. Power density spectra of the normalized lightcurves shown in Figure O In each panel, five power spectra are plotted, one each 
for Models A-E. The panels, from left to right, are the power spectra of the lightcurves in the radio, 1.3 mm, IR, and X-rays. The color 
scheme matches the one we used in Figure [U with the exception that, for the IR and X-rays (the third and fourth panels), we plot the 
purple curve s in da sh as they overlap with the blue curves. In the third panel, the shape of the power spectrum inferred observationally by 
IMever et al ] ll200l) is shown as a dashed black curve, with its uncertainty depicted by a shaded gray band. In the remaining three panels, 
the dashed lines show power-law dependencies that are similar to those exhibited by the simulation data. 
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Figure 8. Cross correlations between the lightcurves at different wavelengths and the lightcurve at 1.3 mm, for the five models shown in 
Figure O For the last panel (labeled as inner X-rays), we consider only the emission that originates from the inner 16M x 16M portion of 
each image. 
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Figure 9. The 1.3mm (black) and IR (green) lightcurves for Model B (neither rescaled nor normalized). The middle and lower rows 
show the 20M x 20M images at 1.3 mm and in the IR at different points in time that correspond to a quiescent reference and to several 
flares. All images use the same quasi-logarithmic color scale (see adjacent color bars), where unity represents the minimum of the peak 
brightness value of the images in each row. The small-amplitude flare in the second column is generated by strong-held gravitational lensing 
of a small, magnetically dominated region that crossed a caustic behind the black hole and gave rise to a nearly perfect Einstein ring. 
The large-amplitude hares in the third and fourth columns are generated by sheared magnetically dominated regions, the light from which 
appears both above and below the equatorial plane, because of gravitational lensing. 
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Figure 10. The radio (3 cm) lightcurve for Model B (neither rescaled nor normalized) and the 256M x 256M images at different points in 
time corresponding to the global minimum of the lightcurve, the two flare events shown in Figure [9] and the global maximum. All images 
use the same quasi-logarithmic color scale (see adjacent color bar), where unity represents the minimum of the peak brightness value of all 
images. The emission originates primarily in the hot, magnetically dominated funnel and is partially obscured by the colder accretion ffow. 
Any gravitationally lensed structures, such as those shown in Figure [9]for the 1.3 mm and IR images, are completely obscured by the cold 
accretion ffow. 


the interplay between emission from a bright funnel and 
osculation by a cold accretion disk. 

4.1.1. Strong-Field Lensing of Magnetic Flux Tubes 

In Figure [9] we plot the 1.3 mm and IR light curves of 
Model B as well as the corresponding images at a qui¬ 
escent reference and during several representative flares, 
with different amplitudes. It is clear from the achromatic 
character of the bright structures in the images as well as 
from their characteristic shapes that strong-field gravita¬ 
tional lensing plays the dominant role in generating these 
flares. 

In some instances, such as the one shown in the sec¬ 
ond column from the left in Figure [H a hot, magneti¬ 
cally dominated structure that has a small spatial extent 
crosses a caustic, giving rise to a nearly perfect Einstein 
ring and an accompanying flux incre ase. This is the phe¬ 
nomen on that was investigated by iRauch fc BlandfordI 
dmi), who calculated the properties of caustics around 
rotating black holes. In other instances, the lensing of 
a sheared flux tube of small radial but larger azimuthal 
extent produces images of this structure to appear both 
above and below the disk and causes large amplitude 
flares (see the right two columns of Figure [9]). This 
phenomenon is very similar to the large increase of the 
brightness of stars during microlensing events. 

Because of their gravitational origin, the presence of 
the large flux excursions in our simulations depends very 
weakly on the simplifying assumptions that are inher¬ 
ent to our model. The only requirement is the fact that 
the flow needs to be optically thin in order for the lensed 
images not to be obscured by intervening material. How¬ 
ever, the actual spectra of the flares depend strongly on 
our prescription for the plasma thermodynamics in the 


magnetically dominated regions. We recall that, in both 
Models A and B, when the plasma-/3 parameter exceeds a 
threshold (here (3 > 0.2) we use a thermal plasma model 
with a constant electron temperature (see discussion in 

IChan et, aUlaTr^ . 

The plasmas in these short-lived, magnetically domi¬ 
nated filaments may not have time to thermalize com¬ 
pletely and this will affect the particular spectrum of 
the flares, predominantly in the IR. Nevertheless, within 
our simplifying assumptions, the IR luminosity at 2.1 fim 
becomes an order of magnitude larger than 1.3 mm lumi¬ 
nosity during the strong flares, from which we can infer 
the spectral index between these two wavelengths as 

^ ^ log(-^2.lMm/i^l.3 mm) ^ Q ^ 

log(2.1 /im/1.3 mm) 

Even though this is not directly comparable to the 
a —0.6 T (k2_value inferred observationally by 
iHornstein et al.l (j2007[ ) for the index of the IR spectrum 
alone, it is significant that our simple thermal model for 
the electrons in the magnetic filaments generate a flare 
spectrum that is substantially bluer than that of the av¬ 
erage emission. 

4.1.2. Bright Radio Funnels Obscured by 
Optically Thick Disks 

In section [331 we asserted that, in the disk-dominated 
SANE models, the large (> 10 hr) positive time lags ob¬ 
tained between the radio and the 1.3 mm light curves are 
artifacts of matching physically unrelated flux variations 
between these wavelengths. Figure [10] shows the radio 
(3 cm) images calculated during the global minimum of 
the lightcurve, the two flare events shown in Figure O 
and the global maximum of the lightcurve. A comparison 
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of Figures [9] and [To] shows that the emission at the two 
wavelengths originates from two different, and discon¬ 
nected locations in the flow. The radio emission comes 
from a hot funnel, which is partially blocked by a colder, 
optically and geometrically thick accretion disk. Note 
also that the scale of these images is much larger than 
the scales of the 1.3 mm and the IR images shown in 
Figure in In contrast, all gravitationally lensed magnetic 
filaments that are responsible for the large 1.3mm/IR 
flares lie very close to the horizon and are completely ob¬ 
scured by the cold accretion disk at radio wavelengths. 

4.2. Jet-dominated MAD models 

The lightcurves of the jet-dominated, MAD models do 
not show any evidence for short-lived flaring events but 
rather slower, correlated flux variations across the radio- 
to-IR part of the spectrum. Moreover, there is marginal 
evidence for the IR light curve leading the 1.3 mm one by 
^ 1 hr and the 1.3 mm lightcurve leading the radio (3 cm) 
one by a comparable amount. 

In this section, we attribute a set of quasi-regular dim¬ 
ming episodes in the 1.3 mm lightcurve of Model C to a 
time-dependent loading of the footpoint of the jet funnel. 
We also point out some evidence for variability generated 
by dissipation of magnetic energy in localized events. 

4.2.1. Time-variable Loading of the Jet Funnel 

The non-spinning MAD Model C shows a peculiar set 
of short timescale, quasi-regular dimming episodes at 
1.3 mm, marked by arrows in the lightcurves plotted in 
Figure Ull The same figure also shows representative im¬ 
ages during the flux minima and recovering phases of two 
such dimming events. It is evident that during each dim¬ 
ming event, e.g., the images shown in the first and third 
panels from left, the bright funnel region becomes less 
active and more collimated. As the brightness recovers, 
the funnel becomes active again, with waves of hotter 
regions propagating along the funnel walls. Note that, 
this effect is quite diffe rent from the radio / mm d imming 
of Sgr A* discussed in lYusef-Zadeh et ah! (j2QlQ[ ). which 
may be caused by obscuration of optically thick plasma. 

Among the five best-fit models of Sgr A* that we are 
considering in this paper. Model C is the only one for 
which the emission is do minated by the fo otpoint of the 
jet funnel (see images in iChan et al.ll2Ql^ . It is, there¬ 
fore, not surprising that this is the only model that is 
characterized by this dimming behavior. The individ¬ 
ual episodes are short-lived (< 1 hr) and quasi-regular 
(once every ^ 1-3hr). Even though the overall behav¬ 
ior of these events is well captured by our simulations, 
the dimming depth is affected by the fact that we excise 
the innermost grid cells near the rotational pole of the 
sp acetime to avoid a rtificially hot pixels, as documented 
in iChan et al ] (1201^ . 

4.2.2. Magnetic Reconnection 

In the rapidly spinning, jet-dominated MAD Models D 
and E, the relatively large extent of the emission region 
across the long-wavelength spectrum and the lack of well 
defined, short-lived flaring events makes it difficult to 
identify the physical origin of the variability in the sim¬ 
ulated lightcurves. Nevertheless, there is some evidence 
for variability generated by dissipation of magnetic en¬ 
ergy in localized events, as shown in Eigure[T2l 


Visual inspection of the images during and after one 
of these events reveals that the increase in the emission 
starts as a localized hot spot in the X-rays, most proba¬ 
bly due to magnetic reconnection and dissipation of the 
magnetic energy. This hot emitting region then spirals 
upwards, and the change in its physical location is mir¬ 
rored in the evolution of its temperature and, hence, in 
the peak wavelength of its emission. This is the origin 
of the ^ 1 hr long time lags between the IR, 1.3 mm, 
and radio lightcurves discussed in section 13.51 This is 
similar to the corona l mass ejecti on episodes observed 
in the Sun (see, e.g., lEorbesII2QQQL for a review) and is 
the closest analogue within G RMHD simulations t o the 
expanding plasmon model of Ivan der LaanI (j 19661 ) that 
has been used in several phenou ienological description of 
time lags fr om Sgr A* (see, e.g.. lYusef-Zadeh et al.ll2QQ8l 
[2001 [ 2 (^ . 

5. DISCUSSION AND CONCLUSIONS 

In this paper, we performed an extensive exploration 
of the predicted variability properties of GRMHD nu¬ 
merical models of Sgr A*, with parameters that have 
been fixed in earlier work (jChan et al.ll201^ by fitting 
their time-averaged properties to the observed spectrum 
of Sgr A* and to the size of its 1.3 mm image. In par¬ 
ticular, we investigated the variability properties of five 
models, two of which were disk-dominated SANE models 
and three were jet-dominated MAD models. 

We found that the disk-dominated SANE models gen¬ 
erated photon-frequency dependent amplitudes of vari¬ 
ability and IR power spectra that are in agreement 
with those inferred observationally. Their lightcurves 
also show a number of short-lived, large flaring events 
that have short (< 1 hr) time lags between the IR and 
1.3 mm, but have no X-ray counterparts. The distribu¬ 
tion of IR flare amplitudes is well described by a log¬ 
normal distribution, with potential excess at large am¬ 
plitudes, in agreement with the observations. We showed 
that strong-field gravitational leasing of hot, magnet¬ 
ically dominated regions is responsible for these large, 
achromatic flares. The reason that such flares were ab¬ 
sent in our earlier work (jChan et al.ll2QQ9[) is because we 
did not include there the effects of gravitational leas¬ 
ing and incorporated a more simplified model for plasma 
thermodynamics. However, we also confirmed our earlier 
result that MHD turbulence alone with thermal plasma 
processes does not produce large X-ray flares, since the 
X-ray emission originates at large radii, where gravita¬ 
tional leasing is negligible. 

We speculate that particle acceleration processes need 
to be incorporated in the plasma model in the GRMHD 
simulations to account for the observed X-ray flaring. We 
expect that the presence of a short-lived distribution of 
relativistic electrons in the inner accretion flow will gen¬ 
erate significant flux excursions from the X-rays down to 
the IR. However, incorporating such non-thermal elec¬ 
trons in the model will not affect significantly the proper¬ 
ties of the flares we report here. Indeed, because the large 
IR flares are the result of strong-held gravitational leas¬ 
ing of very localized magnetic hlaments, the total scatter¬ 
ing optical depth (for Compton emission) or emission vol¬ 
ume (for synchrotron emission) of non-thermal electrons 
in these hlaments will be too small to generate signih- 
cant excursions in the X-ray hux (see also the estimates 
























12 


Time (hr) 



40 20 0 -20 -40 40 20 0 -20 -40 40 20 0 -20 -40 40 20 0 -20 -40 

Relative right ascension (|aas) Relative right ascension (|aas) Relative right ascension (|aas) Relative right ascension (|aas) 


Figure 11. The 1.3mm lightcurve for Model C and the corresponding 20M x 20M images during representative short timescale, quasi¬ 
regular dimming episodes. Because the variation in the flux during each event is small, we use a linear color scale for images (see adjacent 
color bar, where unity represents the minimum of the peak brightness values of the images). Each dimming event occurs when the mass 
loading of the footpoint of the jet funnel is reduced and the funnel itself becomes narrower. 


in lDodds-Eden et~ani2QQ9L section 6). In contrast to the 
SANE models, the variability of the MAD models shows 
only slow trends, with no evidence for short-lived flares at 
any of the wavelengths. The non-spinning jet-dominated 
model, where the 1.3 mm emission is localized at the foot- 
point of the funnel, shows quasi-regular dimming events 
that arise from the time-dependent loading of the funnel. 
The rapidly spinning jet models show weak evidence for 
time lags between wavelengths, with longer wavelengths 
lagging the shorter ones. These time lags seem to arise 
from the upward motion of hot regions heated by mag¬ 
netic reconnection/dissipation events. Overall, the vari¬ 
ability properties of the MAD simulations do not appear 
to agree with those inferred observationally for Sgr A*. 

Our results indicate that, in the near future, it will 
be possible to observe strong-held gravitational lensing 
effects around black holes with the Event Horizon Tele¬ 
scope at 1.3 mm and with GRAVITY in the IR. Even 
though the hare events are too short-lived to obtain 
full images during a single lensing event, their pres¬ 
ence can be inferred from the evolution of the closure- 
phase s obtained along triangles with appropria te base¬ 
lines (jPoeleman et al .112 00 91: 1 Vincent et al.ll2Q14l ). 
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Figure 12. The radio (red)^ 1.3 mm (black), and X-ray (blue) lightcurves for Model E and the corresponding 64M x 64M images after 
a representative reconnection event. The fluxes in the total and inner X-ray lightcurves have been scaled by factors of 10 and 1000, 
respectively, for clarity. The color scales of the images is quasi-logarithmic (see adjacent color bars, where unity represents the minimum 
of the peak brightness values of the images for each row). The first column shows a snapshot where the inner X-ray lightcurve reaches a 
local maximum. We attribute this to a reconnection event, where magnetic energy is dissipated into thermal energy and the topology of 
the held changes. This generates a hot spot in the X-ray image, at the left side of the black-hole shadow. In the second column, a twisted 
feature develops in the left hand side of the 1.3 mm image and moves upwards in the third column, at which point the 1.3 mm lightcurve 
reaches a local maximum. Finally, the fourth column shows a snapshot where the inner X-ray lightcurve reaches a local minimum, but the 
first twisted feature appears in the funnel in the radio image above the cold, obscuring accretion disk. This set of events resembles coronal 
mass ejection episodes observed in the sun. 
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